Figure 5 and 20: 1D plume evolution#
For higher resolution figure, change npos in load-calc-hires-zoom-av.ipynb
%run import-modules-grid
import holoviews as hv
from holoviews import opts
hv.extension('bokeh')
from bokeh.plotting import figure, show
from bokeh.io import output_notebook
!jupyter --version
import os.path
import copy
output_notebook()
Platform: Darwin Kernel Version 24.1.0: Thu Oct 10 21:00:32 PDT 2024; root:xnu-11215.41.3~2/RELEASE_ARM64_T6030
python version: 3.11.10
matplotlib version: 3.9.2
hvplot version: 0.11.0
numpy version: 2.1.2
pandas version: 2.2.3
pickle version: 4.0
yaml version: 6.0.2
pint version: 0.24.3
pyko version: v0.8.3-dev-2024-05-12
print eos_table version: v1.1.5b
Number of CPUs in the system: 12
Selected Jupyter core packages...
IPython : 8.28.0
ipykernel : 6.29.5
ipywidgets : 8.1.5
jupyter_client : 8.6.3
jupyter_core : 5.7.2
jupyter_server : 2.14.2
jupyterlab : 4.2.5
nbclient : 0.10.0
nbconvert : 7.16.4
nbformat : 5.10.4
notebook : not installed
qtconsole : not installed
traitlets : 5.14.3
Read in example vapor plume calculation#
rplumeinitarr = np.asarray([25.e3]) # m
rnebinitarr = np.asarray([250.e6]) #m
tsavearr = np.asarray([1.]) #s
pinitarr = np.flip(np.asarray([80.e9]))
velinitarr = np.asarray([0.]) # m/s
tstall = np.zeros((len(rplumeinitarr),len(pinitarr),len(velinitarr)))
if 0:
fin='./data/temp/vp-h2o-tempgridlong.yml'
fout='./data/temp/vp-h2o-tempgridlong-'
ftype='YAML'
verbose=True
userdtstart=0.0001
usertstepscale=1
binoutput=True
debug=False
initarr=False
#run = RunClass(fin=fin,fout=fout,ftype=ftype) # initialize run parameters class
#
# read in the run template
#readinput_yaml(run,verbose=0)
#print(run.ieos[0].path)
#hugarr = np.loadtxt(run.ieos[0].path+'NEW-SESAME-HUG-E.TXT',skiprows=3,delimiter=',')
#print('Hugoniot file size = ',hugarr.shape)
# vary starting parameters
for ipp in range(len(pinitarr)):
for irp in range(len(rplumeinitarr)):
for ivel in range(len(velinitarr)):
fileid = 'P'+str(np.round(pinitarr[ipp]/1.e9))+'-R'+str(np.round(rplumeinitarr[irp]/1.e3))+'-V'+str(np.round(velinitarr[ivel]/1.e3))
outputfilename = fout+fileid+'.dat'
print("#outputfilename='"+outputfilename+"'")
if os.path.isfile(outputfilename):
%run load-calc-hires-zoom-av.ipynb
print('loading '+outputfilename)
else:
print('file not found: '+outputfilename)
pkopostrack_orig = copy.deepcopy(pkopostrack)
# av version is to mitigate the artificial viscosity error at the material interface in the plot
# note: material interface temperature error suppressed for plotting for +- 5 cells
pkopostrack = copy.deepcopy(pkopostrack_av)
outputfilename='./data/fig5-data.pkl'
with open(outputfilename,"wb") as f:
pickle.dump([pko,pkopostrack,rstall],f)
if 1:
outputfilename='./data/fig5-data.pkl'
with open(outputfilename,"rb") as f:
[pko,pkopostrack,rstall] = pickle.load(f)
jmaxmat1=max(pko['j'][pko['time']==0][pko['mat']==1])
## this tracks the plume front versus time by following the index j at the edge
timeall = np.unique(pko['time'])
imat1 = np.where((pko['mat'] == 1))
jmaxmat1=max(pko['j'][pko['time']==0][pko['mat']==1])
jmat1all =np.where((pko['j']==jmaxmat1))[0]
posmat1all = (pko['pos'][jmat1all])
#print(jmat1all)
#plt.plot(timeall,posmat1all)
#print(timeall,posmat1all)
#edgeline = hv.Curve((timeall,posmat1all)).opts(color='white')
edgelinewtdashthin = hv.Curve((timeall,posmat1all)).opts(color='white',line_dash='dashed',line_width=.5)
edgelinewtdashthick = hv.Curve((timeall,posmat1all)).opts(color='white',line_dash='dashed',line_width=4)
#edgelinebl = hv.Curve((timeall,posmat1all)).opts(color='black')
edgelinebldashthick = hv.Curve((timeall,posmat1all)).opts(color='darkgrey',line_dash='dashed',line_width=4)
edgelinebldashthin = hv.Curve((timeall,posmat1all)).opts(color='darkgrey',line_dash='dashed',line_width=.5)
edgelinebldashthick
#edgelinebldash = hv.Curve((timeall,posmat1all)).opts(color='black',line_dash='dashed',line_width=.5)
print(max(posmat1all))
itmp = np.where(posmat1all == max(posmat1all))[0]
print(itmp)
print(timeall[itmp])
stall_t = timeall[itmp]
stall_p = max(posmat1all)
stall_pt = hv.Scatter((stall_t,stall_p)).opts(size=8, fill_color="grey", line_color="grey")
86.3635503798735
[428]
[11.87500197]
posline = 40*np.power(timeall,1)
possedov = 40*np.power(timeall,.4)
sedovline = hv.Curve((timeall,possedov)).opts(color='darkgrey',line_width=4)
sedovline2 = hv.Curve((timeall+.75,possedov)).opts(color='black',line_width=4)
startline = hv.Curve((timeall,posline)).opts(color='darkgrey',line_width=4,line_dash='dotted')
#edgelinebldashthick*sedovline*startline.opts(ylim=(0,100))
sedovplot=edgelinebldashthick*sedovline*sedovline2*startline*stall_pt.opts(ylim=(0,100),width=500,height=400,xlabel='Time (hr)',ylabel='Radial Position (Mm)',fontsize={'title': 16, 'labels': 14, 'xticks': 12, 'yticks': 12, 'cticks':12})
#sedovplot=edgelinebldashthick*sedovline*sedovline2*startline*stall_pt.opts(ylim=(0,100),width=1200,height=800,xlabel='Time (hr)',ylabel='Radial Position (Mm)',fontsize={'title': 26, 'labels': 24, 'xticks': 22, 'yticks': 22, 'cticks':22})
sedovplot
# Figure 20 was saved as a png from the widget below
hv.save(sedovplot,'./plots/Fig20-sedovplot.png')
# for each initial position, track
# 0 - time of first nebula shock
# 1 - time of first plume front
# 2 - peak velocity nebula shock and time
# 3 - peak velocity of plume shock and time
# 4 - peak temperature of nebula shock and time (not the same as velocity)
# 5 - peak temperature of plume shock and time (maybe same as velocity)
# 6 - peak pressure of nebula shock and time (not the same as velocity)
# 7 - peak pressure of plume shock and time (maybe same as velocity)
# 8 - time with first shock
npos = 1024
waveinfoarr = np.zeros([9,npos])
waveinfoarr[6,:]=np.log10(1.4064186900465516E-07)
tind=[]
tind2=[]
tracktimeall = np.unique(pkopostrack['time'])
ntime = len(tracktimeall)
itimestall = np.where(tracktimeall > tstall/3600.)[0][0]
trackposall = np.unique(pkopostrack['pos'])
iposstall = np.where(trackposall < rstall/1.e6)[0]
for ipos in np.arange(2,len(trackposall)):
tind=[]
tind2=[]
# first motion?
pkoposone = pkopostrack[pkopostrack['pos']==trackposall[ipos]]
pkoposonemat1 = pkoposone[pkoposone['mat']<1.5]
pkoposonemat2 = pkoposone[pkoposone['mat']>1.5]
ind1= np.where(np.asarray(pkoposone['up']>10))[0] # moving positions
if len(ind1)>0:
waveinfoarr[0,ipos]=tracktimeall[ind1[0]]
#waveinfoarr[8,ipos]=
#print(tracktimeall[ind1[0]])
ind2 = np.where(np.asarray(pkoposone['mat']==1))[0] # mat1 positions
if len(ind2)>0:
waveinfoarr[1,ipos]=tracktimeall[ind2[0]]
#print(ind1[0],ind2[0])
if len(ind1)>0:
tind = np.arange(ind1[0],ind2[0]-1)
tind2 = np.arange(ind2[0],itimestall)
else:
if len(ind1)>0:
tind = np.arange(ntime)
else:
tind=[]#np.arange(npos)
if len(tind)>0:
#waveinfoarr[2,ipos] = np.max(np.asarray(pkoposone['up'])[tind])
#waveinfoarr[4,ipos] = np.max(np.asarray(pkoposone['temp'])[tind])
waveinfoarr[2,ipos] = np.max(pkoposonemat2['up'])
waveinfoarr[4,ipos] = np.max(pkoposonemat2['temp'])
waveinfoarr[6,ipos] = np.max(pkoposonemat2['rho'])
if len(tind2)>0:
#waveinfoarr[3,ipos] = np.max(np.asarray(pkoposone['up'])[tind2])
#waveinfoarr[5,ipos] = np.max(np.asarray(pkoposone['temp'])[tind2])
waveinfoarr[3,ipos] = np.max(pkoposonemat1['up'])
waveinfoarr[5,ipos] = np.max(pkoposonemat1['temp'])
waveinfoarr[7,ipos] = np.max(pkoposonemat1['rho'])
s1len = len(np.where(timeall<17)[0])
s1posall = np.zeros([s1len])
for it in np.arange(s1len):
temp = pkopostrack[pkopostrack['time']==timeall[it]]
temp2 = temp[temp['temp']>250]
#print(it,min(temp2['pos']))
s1posall[it]=(min(temp2['pos']))
# display(temp2)
# print(temp[pkopostrack['temp']>400][0])
#plt.plot(timeall[0:s1len],s1posall)
s2frontblthick = hv.Curve((timeall[15:s1len],s1posall[15::])).opts(color='black',line_dash='dotted',line_width=4)
s2frontwtthick = hv.Curve((timeall[15:s1len],s1posall[15::])).opts(color='black',line_dash='dotted',line_width=4)
s2frontblthick
#s1front = hv.Curve((timeall,waveinfoarr[2,:])
print(len(s1posall))
print(len(timeall[0:s1len]))
505
505
s1len = len(np.where(timeall<6.3)[0])
s1posall = np.zeros([s1len])
for it in range(s1len):
temp = pkopostrack[pkopostrack['time']==timeall[it]]
temp2 = temp[temp['up']>0.8]['pos']
if len(temp2)>1:
s1posall[it]=max(temp2)
# display(temp2)
# print(temp[pkopostrack['temp']>400][0])
#plt.plot(timeall[0:s1len],s1posall)
s1frontblthick = hv.Curve((timeall[0:s1len],s1posall)).opts(color='black',line_dash='solid',line_width=4)
s1frontwtthick = hv.Curve((timeall[0:s1len],s1posall)).opts(color='white',line_dash='solid',line_width=4)
s1frontblthick
#s1front = hv.Curve((timeall,waveinfoarr[2,:])
## this tracks the plume front versus time by following the index j at the edge
jradstall2=min(pko['j'][pko['time']==0][pko['pos']>rstall/1.e6*1.25])
print(jradstall2)
streamline0 = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='darkgrey',line_dash='solid',line_width=1.5)
streamline0wt = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='white',line_dash='solid',line_width=1.)
jradstall2=min(pko['j'][pko['time']==0][pko['pos']>rstall/1.e6])
print(jradstall2)
streamline1 = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='darkgrey',line_dash='solid',line_width=1.5)
streamline1wt = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='white',line_dash='solid',line_width=1.)
jradstall2=min(pko['j'][pko['time']==0][pko['pos']>rstall/1.e6*.75])
print(jradstall2)
streamline2 = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='darkgrey',line_dash='solid',line_width=1.5)
streamline2wt = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='white',line_dash='solid',line_width=1.)
jradstall2=min(pko['j'][pko['time']==0][pko['pos']>rstall/1.e6*.5])
print(jradstall2)
streamline3 = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='darkgrey',line_dash='solid',line_width=1.5)
streamline3wt = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='white',line_dash='solid',line_width=1.)
edgelinebldashthick*streamline1*streamline2*streamline3*streamline0
543
475
407
337
# bigger fonts
# extract data from pko dataframe for an x-t diagram
fs1=18 # title
fs2=18 # labels
fs3=18 #xticks,yticks
fs4=16 #cticks
xrr=(0,14)
yrr=(0,100)
wpix = 800
hpix = 500
heatmap1 = pkopostrack.hvplot.heatmap(y='pos', x='time', C='up',
title="A. Velocity Evolution",
clabel='Radial velocity (km/s)',clim=(-4,4),
xlabel='Time (hr)',xlim=xrr,ylim=yrr,
ylabel='Radial Position (Mm)',
width=wpix, height=hpix).opts(cmap='bwy',fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4})
heatmap1b = pkopostrack.hvplot.heatmap(y='pos', x='time', C='up',
title="A. Vapor Plume Features",
clabel='Radial velocity (km/s)',clim=(-4,4),
xlabel='Time (hr)',xlim=xrr,ylim=yrr,
ylabel='Radial Position (Mm)',
width=wpix, height=hpix).opts(cmap='bwy', fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4})
heatmap2 = pkopostrack.hvplot.heatmap(y='pos', x='time', C='pres',
title="B. Pressure Evolution",
clabel='log₁₀P (Pa)',
xlabel='Time (hr)',
ylabel='Radial Position (Mm)',xlim=xrr,ylim=yrr,
clim=(-5,1),
width=wpix, height=hpix).opts(cmap='bgy', fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4})
heatmap2b = pkopostrack.hvplot.heatmap(y='pos', x='time', C='pres',
title="B. Distinct Zones Followed by Mixing",
clabel='log₁₀P (Pa)',
xlabel='Time (hr)',
ylabel='Radial Position (Mm)',xlim=xrr,ylim=yrr,
clim=(-5,1),
width=wpix, height=hpix).opts(cmap='bgy', fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4})
heatmap3 = pkopostrack.hvplot.heatmap(y='pos', x='time', C='rho',
title="C. Density Evolution",
clabel='log₁₀ρ (g/cm³)',
xlabel='Time (hr)',
ylabel='Radial Position (Mm)',xlim=xrr,ylim=yrr,
width=wpix, height=hpix,clim=(-9,-5)).opts(cmap='bgy', fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4})
heatmap3b = pkopostrack.hvplot.heatmap(y='pos', x='time', C='rho',
title="C. Size Sorting",
clabel='log₁₀ρ (g/cm³)',
xlabel='Time (hr)',
ylabel='Radial Position (Mm)',xlim=xrr,ylim=yrr,
width=wpix, height=hpix,clim=(-9,-5)).opts(cmap='bgy', fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4},xlim=(0,13))
heatmap4 = pkopostrack.hvplot.heatmap(y='pos', x='time', C='temp',
title="D. Temperature Evolution",
clabel='T (K)',
xlabel='Time (hr)',
ylabel='Radial Position (Mm)',xlim=xrr,ylim=yrr,
width=wpix, height=hpix,clim=(200,2000)).opts(cmap='plasma', fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4})
heatmap4b = pkopostrack.hvplot.heatmap(y='pos', x='time', C='temp',
title="D. Thermal and Redox Zones",
clabel='T (K)',
xlabel='Time (hr)',
ylabel='Radial Position (Mm)',xlim=xrr,ylim=yrr,
width=wpix, height=hpix,clim=(200,2000)).opts(cmap='plasma', fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4})
linegroup=edgelinebldashthick*stall_pt*s2frontblthick*s1frontblthick
linegroupwt=edgelinewtdashthick*stall_pt*s2frontwtthick*s1frontblthick
streamlines=streamline1*streamline2*streamline3*streamline0
streamlineswt=streamline1wt*streamline2wt*streamline3wt*streamline0wt
#heatmap*edgeline
#group=(heatmap*edgelinedash+heatmap1*edgelinebldash+heatmap2*edgelinedash+heatmap3*edgelinedash).cols(2)
grouplabels2=(heatmap1b*linegroup*streamlines+heatmap2b*linegroupwt*streamlineswt+heatmap3b*linegroupwt*streamlineswt + heatmap4b*linegroupwt*streamlineswt).cols(2)
grouplabels2